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We present experimental results on a model system for studying wave propagation in a complex 
medium exhibiting low frequency resonances. These experiments enable us to investigate a funda- 
mental question that is relevant for many materials, such as metamaterials, where low-frequency 
scattering resonances strongly influence the effective medium properties. This question concerns the 
effect of correlations in the positions of the scatterers on the coupling between their resonances, and 
hence on wave transport through the medium. To examine this question experimentally, we measure 
the effective medium wave number of acoustic waves in a sample made of bubbles embedded in an 
elastic matrix over a frequency range that includes the resonance frequency of the bubbles. The 
effective medium is highly dispersive, showing peaks in the attenuation and the phase velocity as 
functions of the frequency, which cannot be accurately described using the Independent Scattering 
Approximation (ISA). This discrepancy may be explained by the effects of the positional correla- 
tions of the scatterers, which we show to be dependent on the size of the scatterers. We propose a 
self-consistent approach for taking this "polydisperse correlation" into account and show that our 
model better describes the experimental results than the ISA. 

PACS numbers: 43.35.+d, 43.20.+g 



I. INTRODUCTION 

Propagation of waves in complex media continues to be 
a very active subject of research. Of particular interest 
are complex media containing scatterers with resonances 
at low frequencies. Indeed, when the wavelength is large 
compared to the typical size of the inhomogcncities, one 
can use an effective medium approach, i.e., consider that 
the wave propagates as it would do in a homogeneous 
medium with an effective wavevector whose value is re- 
lated to the composition and the structure of the mate- 
rial. In most cases, the effective medium has properties 
that are close to the average of that of its components. 
But when there are resonant scatterers, their contribu- 
tion to the wave field can be strong enough to signifi- 
cantly alter the propagation. In some cases, the effective 
media show intriguing properties, such as a negative re- 
fractive index, which are not encountered in nature; such 
materials are called metamaterials [lrlij]. 

One of the key issues for wave propagation in strongly 
scattering materials is the effect of coupling between the 
scatterers near resonance, an effect that makes the tradi- 
tional independent scattering approximation (ISA) inad- 
equate. Indeed, when the resonances are strong and/or 
the concentration of scatterers high, positional corre- 



lations and multiple scattering loops may have a non- 
negligible contribution. To gain a better understand- 
ing of the underlying physics, experiments on a well- 
controlled system with strong resonances are needed. 
One such system can be found in acoustics, where bubbly 
liquids are often regarded as model systems for the study 
of wave propagation with low frequency resonances. In- 
deed, air bubbles in a liquid have a particularly low reso- 
nance, known as the Minnaert resonance, with an angular 
frequency given by 



R 



where R is the bubble radius, v\ is the velocity of longi- 
tudinal waves in air, and pi and po are the mass densities 
of the air and the liquid, respectively. As the ratios of 
the densities and the velocities are small, it is easy to see 
that this resonance is at very low frequency: 
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= \/3pi/po— < 1. 

Here, fo is the longitudinal sound velocity in the liquid. 
This equation shows that, at resonance, the wavelength 
is large compared to the bubble size. This is important 
because, then, at the resonance frequency, not only is the 
response of the bubble strong, but several bubbles (and 
potentially many) can be driven in phase, thus yielding a 
collective constructive contribution to the pressure field. 
This explains why a minute quantity of gas bubbles can 
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dramatically change the effective acoustic properties of a 
liquid. 

If the bubble is no longer in a liquid, but in an clastic 
medium with shear velocity uq, its resonance frequency 
becomes @, 0: 
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so the same condition of small radius with respect to 
wavelength persists as long as the elastic medium is a 
soft solid in which «o C wo- As a result, elastic bub- 
bly media are also good candidates as model systems for 
investigating wave propagation in the presence of low fre- 
quency resonances. In addition, they offer two practical 
advantages compared to bubbly liquids: (1) bubbles do 
not move, which means that a precise knowledge of their 
positions can be obtained; (2) the frequency of the reso- 
nance can be shifted by tuning the shear velocity in the 
clastic matrix. Bubbly soft media with negligible shear 
velocity have already been studied [7 9], and the ISA 
was found to reliably predict the experimental results, 
at least for concentrations of bubbles up to 1%. The 
present article focuses on a bubbly polydimethylsiloxane 
(PDMS) sample, in which the shear modulus is expected 
to significantly change the bubble's response to the inci- 
dent ultrasonic field and the propagation of waves in the 
medium. Our aim is to see if the ISA is still relevant in 
this case. Beyond the fundamental aspect of the study, 
and the potential applications for design of metamateri- 
als, it should be noted that a better understanding of the 
acoustic properties of bubbly elastic media is interesting 
per se, since many practical applications, in the medical 
or the food science fields for instance, involve the pres- 
ence of bubbles in a soft elastic matrix. 

The present paper is organized as follows. The next 
section gives a brief overview of the ISA, with a focus on 
its application to acoustic propagation in bubbly media. 
Then we describe how our bubbly samples were made 
and characterized. Section IIVI reports the experimental 
measurements of velocity and attenuation in one typi- 
cal sample, and compares the experimental results with 
the ISA. Two main limitations of the model are also dis- 
cussed: the assumption of sample homogeneity on length 
scales larger than the scatterer size (e.g., gradients in the 
concentration of the scatterers) and the neglect of cor- 
relations in the positions of the scatterers. Methods for 
overcoming these limitations are proposed. 



II. THE INDEPENDENT SCATTERING 
APPROXIMATION 

Based on a diagrammatic representation of scatter- 
ing events the ISA predicts the following effective 
wavenumber k for a medium with n scatterers per unit 
volume defined by their scattering function / s : 



where ItQ is the wavenumber in the pure medium (ma- 
trix). The ISA is actually equivalent to the model de- 
veloped, with another approach, by Foldy in 1945 fTTj ) . 
In the following, we will use equally "ISA" or "Foldy's 
model" when refering to equation ([2]). 

In eq. ©, the effect of the scatterers on propagation 
in the pure medium is described by a term proportional 
to nf s , i.e., a simple addition of the individual scattering 
events, which are assumed to be independent. Models 
have been proposed for incorporating either the correla- 
tions between the scatterings [l2| or the loops, i.e., the 
events in which the same scatterer or spot inside the sam- 
ple is visited more than once by the wave [l3|, [lj[ . We 
discuss only the ISA in this section, since it is the sim- 
plest model and it gives reasonable predictions of wave 
transport in many cases. We will introduce later, in sec- 
tion IIV CI and the appendix, extensions to the ISA for 
dealing with the positional correlations of the scatterers. 

Equation ^ considers a monodisperse assembly of 
scatterers. In practical applications, one is often dealing 
with polydispcrsc media, in which case the ISA predicts 
that 



4irn{R)dRf s (R), 
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where n(R)dR is the number of bubbles per unit volume 
whose radius is between R and R + dR, and f s (R) is the 
scattering function for a scatterer of radius R. 

Let us apply the ISA to a bubbly medium. When the 
scatterer is a bubble, its isotropic scattering function is 
given by 
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where 6 is the damping constant due to thermal, viscous 
and radiative losses [15j and ujm the resonance angular 
frequency given by Eq. ([!]), which can be written 
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with k the polytropic index for the transformations un- 
dergone by the gas [lj| (k = 1.4 for air), Pi the pres- 
sure of the gas in the bubbles, po the mass density of 
the matrix, and fx' the real part of the shear modulus of 
the matrix. Surface tension effects are negligible for the 
bubble size considered in our experiments. Mode con- 
version is not considered in this model since it has been 
shown to be negligible for this system on account of the 
large difference between the shear modulus of the ma- 
trix and the longitudinal one (uq <C vo) [lj|. Hence the 
only significant contribution of the nonzero shear modu- 
lus is its effect on the resonance frequency of the bubble, 
which is increased by the stiffness of the matrix. Another 
approximation we make is to reduce the scattering func- 
tion to its isotropic component, neglecting the dipolar 
and higher contributions. This approximation is gener- 
ally good when the wavelength is large compared to the 
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size of the scatterers, and it will be shown to be accurate 
for describing our experiments. 

In figure [U we give two examples of the effective at- 
tenuation and velocity as functions of the frequency, as 
calculated by Foldy's model. For the sake of generality, 
we consider log-normal distributions of bubble size: 
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where Rq is the median radius, e the width of the distri- 
bution, and no the bubble concentration. Note that for 
distributions that are not too broad (e < 0.2), log-normal 
distributions are close to normal distributions: then Rq 
and e can be viewed as the mean radius and coefficient of 
variation (i.e., standard deviation divided by the mean), 
respectively. The latter quantity, which we will refer to 
as the polydispersity factor, measures the range of bub- 
ble sizes in the distribution. The bubble concentration 
no is related to the gas volume fraction $ by 
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The bubble size distribution in Fig. [T] is described by 
Ro = 150 urn, e = 25% and $ = 2%. Two host me- 
dia are considered: water (dashed lines), and an elastic 
matrix with u' = 0.7 MPa, a" = 0.4 MPa (solid lines). 
For both media, very dispersive behavior is predicted: 
the attenuation and velocity show peaks, which occur at 
higher frequencies in the elastic medium. 

At low frequencies (oj < ojm), sound propagation in 
the bubbly medium is weakly attenuated, and its very 
low phase velocity is accurately predicted by Wood's re- 
lation [l7|, modified to account for the elastic effects. For 
5xl0 -4 < $ < 10 _1 , this relation can be approximated 
(within 5 %) as 
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where P\ is the pressure of the gas in the bubbles and po 
the density of the matrix. Eq. ((5J) gives, in our example, 
Vu = 122 and 248 m/s in water and the elastic matrix, 
respectively. 

Near the Minnaert frequencies of the bubbles, the 
acoustic waves are highly attenuated by the bubbly 
medium: the attenuation reaches a peak a max at fre- 
quency /max, corresponding to the resonance frequency 
for a monodisperse distribution. In water, the peak is 
very sharp and associated with the bubble resonances 
(/max = 22 kHz as expected by Eq. (|5|) for a 150-/im- 
radius bubble). Eqs. © and (dJ) give the physical origin 
of this sharp peak: at resonance, the scattering function 
f s is large and imaginary, which leads to a large imag- 
inary part of the effective wave number k. This sharp 
resonant effect is followed by a broader regime (roughly 
60-400 kHz in our example) in which the attenuation is 
still large. This is explained by the negative response of 
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FIG. 1: (Color online) Predictions of the ISA for the atten- 
uation a (top) and phase velocity v (bottom) as functions of 
frequency for bubbles in water (dashed lines) and in an elas- 
tic matrix (solid lines) with \i! = 0.7 MPa, fi" — OA MPa; the 
bubble concentration is 2 % and the mean radius is taken to be 
150 fim with a 25 % polydispersity. The characteristic features 
(fmax, Oi m ax, c*hf, vjf, and v hj ) are indicated. 



the bubbles. Indeed, when driven at a frequency higher 
than its resonance frequency, a bubble has a scattering 
function which reduces approximately to f s = —R when 
S <C l(see Eq. (QJ). Even though the bubble is not res- 
onating, its contribution is far from negligible. Indeed, 
one can rewrite Eq. © as 



k* = ki l - 



RX 2 
nd 3 



(9) 



where d is the typical distance between bubbles. As 
long as the wavelength is large compared to this dis- 
tance (X/d 1), the square of the effective wavenum- 
bcr is negative: waves are evanescent. In our example of 
2 % bubbly water with 150 ^m-radius bubbles, the typ- 
ical distance between bubbles is 0.9 mm, which means 
that RX 2 /(nd 3 ) > 1 for / < 380kHz. 

In our example of an elastic medium, the high viscos- 
ity (n" of 0.4 MPa corresponds, at 0. 1 MHz, to a viscosity 
640 times larger than that of water) tends to damp the 
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oscillations of the bubbles at resonance, so the attenu- 
ation peak is lower and less sharp (see figure [1} . As a 
paradoxical consequence, the more viscous the medium 
in which bubbles are embeded, the lower the maximum of 
attenuation of sound. Note that the sharp peak of atten- 
uation also dissapears when the bubble size distribution 
is too broad. Interestingly, even if the resonant effects are 
too weak to give rise to significant resonant attenuation, 
they still lead to an evanescent regime. In this case, the 
position of the attenuation peak does not correspond to 
the resonance of the average size bubbles, but rather to 
a frequency that is just above the resonance frequency of 
the small bubbles of the distribution, as it corresponds to 
all the bubbles being driven at a frequency higher than 
their resonance frequencies. For instance, in our exam- 
ple /max = 90 kHz for the bubbly elastic medium, which 
corresponds to a radius of 100 /im according to Eq. ([5]). 

When the resonances are smoothed out, an approxi- 
mate formula can be established for the maximum of the 
attenuation: 



(10) 



2V3$ 



J R exp(2e 2 



which has the benefit of depending only on the bubble 
concentration and size distribution (not on either the 
damping constant of the bubbles or the shear modulus of 
the matrix). 

At higher frequencies, the evanescent regime disap- 
pears: the attenuation and velocity reach constant val- 
ues, which can be approximated by: 
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(11) 
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Interestingly, if Foldy's model applies, one can deter- 
mine the 5 parameters $, Rq, e, [i! , and vq from the 
measurement of the 5 quantities v\{, /max, oi max , cthf, and 
Vhf in a bubbly medium. 



III. MATERIALS AND METHODS 

A. Sample preparation 

Samples were prepared by mixing 100 g of RTV 615 
(GE Silicones, 90 g monomers and 10 g hardener) in a 
600 mL beaker with a Sorvall Omni Mixer (Ivan Sorvall, 
Inc., Norwalk, CT; nominal operating speed 16,000 rpm) 
for 1 minute to incorporate bubbles. Then, after waiting 
for an interval of 10-20 minutes for the biggest bubbles 
to rise and disappear, the bubbly liquid was poured into 
a cell (thickness 2.6mm, diameter 11cm), which was ro- 
tated at 5 rpm to prevent creaming (bubbles rising to 
the top of the container). At room temperature, a solid 
sample was obtained in 1 day. 




FIG. 2: (Color online) Example of an image taken for the 
bubble size analysis (optical characterization). 



B. Sample characterization 

Testing of Foldy's model requires that some of the pa- 
rameters of the sample are known: the bubble size distri- 
bution (n(R)), the matrix rheology (/i) and the velocity 
of sound in the matrix (vq). 

a. Optical characterization was made by taking pic- 
tures at different positions in the sample. The depth of 
field was larger than the thickness of the sample, which 
means that the total volume of the image was known 
(16 x 11 x 2.6 mm 3 ). A typical image is shown in fig- 
ure[2j the contrast was good enough for an automatic size 
analysis to be performed (ImageJ). The radius measure- 
ment was possible with a one pixel accuracy, i.e., 5/xm. 
A total of 839 bubbles was analyzed. To avoid biased 
measurements of big bubbles, overlaid bubbles were ex- 
cluded from the size analysis. However, the total number 
of bubbles was determined (1098) for a correct estima- 
tion of the bubble volume fraction, assuming the same 
size distribution for the overlaid bubbles. The radius 
distribution was centred around 165 ± 5 /jm with a poly- 
dispersity of 21 ± 1 % and a volume fraction of 2 ± 0.2 % 
(see histogram in figure [3]). 

b. Rheological characterization was necessary to 
measure the shear modulus over the frequency range used 
for the ultrasonic experiments (30-800 kHz). We used a 
shear wave reflection technique (see [HI or [l!| for de- 
tails) to measure the complex shear modulus /i = /i' + i/i" 
from 300 to 500 kHz. For lower frequencies (around 
30 kHz) an estimation of /j, was possible by measuring 
the acoustic response of a single bubble (see [20] for a 
description of the technique). Over the 30-500 kHz fre- 
quency range, our measurements are consistent with the 
following behaviour of the real (//) and imaginary (//") 
parts of the shear modulus: 



/j = 0.6 + 0.7 x /, 
fi" = 0.2+1.8 x /, 



(13a) 
(13b) 



where the moduli are in MPa and the frequency / in 
MHz. Equations ([TU)) only give the order of magnitude 
of the shear modulus for a given RTV615 sample. In- 
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FIG. 3: (Color online) Histogram of the bubble radius dis- 
tribution for the sample in figure [Jl obtained from the image 
analysis. The superposed curve (circles plus connecting lines) 
corresponds to the result of the x-ray tomography analysis. 



deed, the precise protocol used for preparation, such as 
the temperature and time of the hardening process for 
instance, may have a non-negligible impact on the mod- 
ulus plj . As a consequence, deviations from the predic- 
tions of Eqs. (|13p can be observed for a given sample. 

The stationary limit (/ = MHz) of our estimate of 
fi' is in good agreement with the value reported in the 
literature plj . 

c. Ultrasonic characterization was also performed 
by standard techniques 0] to measure the velocity and 
attenuation of sound in the pure PDMS. The longitudi- 
nal velocity was found to be constant over the frequency 
range investigated (30kHz - 5 MHz): vq = 1.02mm//is, 
in accordance with previously reported values [22| . For 
the attenuation, the following law was found: «o = 
0.023/ 154 mm- 1 (with frequency / in MHz). 

d. X-ray characterization. To investigate the spatial 
distribution of the bubbles in one of the samples, x-ray 
tomography images were acquired, using an Xradia Mi- 
croXCT 3D x-ray transmission microscope (18/im/px). 
After the ultrasonic measurements, the sample was cut 
in pieces small enough to fit in the apparatus (7.5 x 10 x 
2.6mm 3 ), and 21 pieces were imaged, for a total of 3700 
bubbles. A software program was written to enable the 
reconstruction of the 3D structure of each piece: for ev- 
ery bubble i, its radius R\ and position [xt, yi, Zi] were 
determined with 1 pixel accuracy. The size distribution 
was determined and successfully compared to the results 
of the optical characterization: we found a mean radius 
of 160 ± 18 /im with a polydispersity of 25 ± 1 % and a 
volume fraction of 2.2 ± 0.3% (see figure [3] and table [I]). 
Moreover, the 3D data gave insight into the homogeneity 
of the sample, as well as on the correlations in the bubble 
positions, as described in sections llV Bl and IIV CI 



TABLE I: Parameters of the bubble size distribution, extracted 
from the optical and x-ray measurements. 





volume fraction 
$ (%) 


average radius 
(R) (urn) 


polydispersity 

(R) ( /o > 


optics 
x-ray 


2.0 ±0.2 
2.2 ±0.3 


165 ±5 
160 ± 18 


21 ± 1 
25 ± 1 



C. Ultrasonic measurements 

The acoustic properties of the samples were measured 
with the following set-up. In a large tank (60 x 60 x 
120 cm 3 ) filled with reverse osmosis water, a piezoelectric 
transducer generated a pulse that propagated through 
water, traversed the sample and was detected by a hy- 
drophone. Because the attenuation in the sample was 
large, and because the divergence of the beam was not 
negligible (especially at low frequencies), the use of a 
screen was essential to reduce spurious signals. The aper- 
ture of the screen (6 cm) was larger than the wavelength 
of the pulse (3.75 cm at 40 kHz) to limit diffraction effects, 
but smaller than the diameter of the sample (11 cm). 

Gaussian pulses, with central frequencies ranging from 
40 to 600 kHz, were generated by an Arbitrary Wave Gen- 
erator, and three different transducers, having central fre- 
quencies of 100, 250 and 500 kHz, were used to cover the 
range of frequencies. The pulses were recorded, with a 
hydrophone, in two different cases: when the sample was 
mounted on the screen (si(i)), and when the sample was 
absent (s2(i)). The signals were averaged over 100 ac- 
quisitions when the attenuation was low (for reference 
measurements, for example), and up to 5000 acquisitions 
for highly attenuated signals. Signals si(t) and S2{t) were 
then truncated to eliminate spurious echoes, and Fourier 
transformed into S\(uj) and S^O^) respectively. Then, 
T(oj), the ratio of the transmission with and without the 
sample in the path of the acoustic beam at a given an- 
gular frequency, was calculated. 

For an incoming plane monochromatic wave exp(ifcir — 
iujt), T(w) is given by 
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where k and k w are the wave numbers in the sample 
and in the water respectively, Z and Z w the acoustic 
impedances, and d the thickness of the sample. Because 
the impedance Z depends on k (Z = ptu/k, with p the 
mass density of the sample), equation (fT4"]) is not directly 
invertible to extract k: an iteration method was used [8|. 
From k, the phase velocity v = w/Re(fc) and the atten- 
uation a = 2Im(fc) can be calculated. Error-bars on the 
measurements were evaluated by taking into account the 
following sources of uncertainty: thickness of the sample, 
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noise in the Fourier transforms, truncation of the time 
pulse. 

An important issue in measuring the acoustic proper- 
ties of a medium is the statistical relevance of results for 
a single sample. Theories predict the average properties, 
i.e., what one can measure when averaging over different 
samples having the same general parameters. To obtain 
insight on this issue, we measured the transmitted signal 
not only in the centre of the acoustic beam, but also on a 
7x7 grid with a 5 mm step. From this scanning, a "typi- 
cally averaged signal" could be obtained (with correction 
for the geometrical differences), and this was compared 
to the centre pulse. 



IV. EXPERIMENTAL RESULTS AND 
DISCUSSION 



Several bubbly PDMS samples were tested, with dif- 
ferent void fractions (0.5-3.5 %) and polydispcrsitics 
(20-50%), but with the same average radius (around 
150 /mi). In this article, we focus on representative re- 
sults obtained with one selected sample, which was char- 
acterized by x-ray tomography Figure 2] shows the 
measured attenuation and phase velocity in the sam- 
ple as functions of the frequency. The general aspect 
of the curves follows what one would expect for a bubbly 
medium (see figure [TJ : both the attenuation and the ve- 
locity exhibit peaks in frequency with very large values at 
their maxima (3.5 mm -1 and 4.5mm//ts, respectively), 
the velocity is low (0.26 mm//is) at low frequencies, and 
close to the value in the pure PDMS at higher frequencies. 
Note that two cases are considered: when only one pulse 
is analysed (open circles) and when the average is made 
over the 49 positions of the hydrophone (solid circles). 
Except for the peak in velocity, both analyses give simi- 
lar results. It indicates that substantial information can 
be obtained from a single pulse measurement, which is in- 
teresting for practical purposes. Note that measurements 
are difficult when the velocity is high because small time 
shifts have to be extracted from the phase shift which 
is dominated by the complex impedance rather than the 
transit time. 

The five typical quantities introduced in section |TT] 
can be measured and compared to equations ([5]), and 
([8])-(fl2j). calculated with the parameters found in sec- 
tion IIII Bl As shown in table HU the agreement is rather 
satisfying: / max and a^f are within 10% and the other 
quantities within 20%. Note that for v\{ and vm the mea- 
sured values may not correspond exactly to the asymp- 
totic values because of the limited range of frequencies in 
the measurements. 



A. Comparison with ISA 

Figurc[5]offers a further step in the comparison process. 
Predictions of the ISA model are plotted as dotted lines. 
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FIG. 4: (Color online) Attenuation a and phase velocity v in 
the bubbly PDMS sample, measured from a single pulse ( open 
circles) and from an average pulse (solid circles). The velocity 
in the pure PDMS is shown by the solid line. The attenuation 
in pure PDMS is, on average, 2 orders of magnitude less than 
in the bubbly sample (too small to be visible on the graph). 



TABLE II: Five characteristic quantities measured experimen- 
tally (first line), and predicted by the model (second line). 





/mas 


Qimax 


Qhf 




VhS 




(kHz) 


(mm" 1 ) 


(mm" 1 ) 


(m/ S ) 


(m/s) 


exp. 


100 


3.5 


0.40 


260 


1180 


mod. 


90" 


2.8 


0.35 


224 


1020 



"Calculated with Eq. ([5} by taking the radius of the small bubbles 
of the distribution (i.e., 100 fivn, see figure [3}. This is appropriate 
since for polydisperse distributions of sizes, the peak in attenuation 
corresponds to the resonances of the smaller bubbles, as noted in 
section Ull 



Parameters for the model were taken as measured in sec- 
tion IIII BL the histo gram measured in the x-ray experi- 
ment being chosen for n(R). In the low and high limits of 
the frequency range, the comparison is good. However, 
the attenuation and velocity peaks are not correctly pre- 
dicted: their frequencies and amplitudes are both smaller 
than the experimental data. Note that this discrepancy 
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was observed for all the samples we investigated. The 
comparison is slightly better when the dipole terms in 
the bubbles' responses are considered (solid lines) but 
the discrepancy in the peaks' positions and magnitude is 
still large. Higher order scattering has negligible contri- 
butions. 
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FIG. 5: (Color online) Close-up of the attenuation and ve- 
locity peaks. Solid circles correspond to experimental data 
(same as in Fig.^. Comparison with Foldy's model is shown 
for monopolar scattering of the bubbles (dotted lines) and 
monopolar + dipolar scattering (solid lines). Note that there 
are no adjustable parameters in this comparison of the the- 
ory with experiment. Results of fitting the model to the data 
by allowing either the size distribution (dashed lines) or the 
shear modulus (triangles) to vary are also shown. In the in- 
set, the size distribution determined by the x-ray experiment 
(solid line) is compared to the one inferred from the fit to the 
ultrasonic data (dashed lines). 



To check whether the discrepancy could be explained 
by the uncertainty in the measured size distribution, we 
determined the n(R) function that was needed to achieve 
a better fit. As shown in the inset of Fig. [SJ this leads 
to a totally unrealistic distribution, with smaller bub- 
bles (mean radius 120 /im), a narrower distribution (10 % 
polydispersity) , and a smaller void fraction (1.1 %). One 
can understand why such a size distribution is obtained 
by looking at equation @. If the peaks are at too low fre- 
quencies and with too low amplitudes, it means that the 
nf s term is not high enough compared to i.e., that the 
bubble contribution is too weak. A way of fixing this is to 
consider smaller bubbles, in order to shift the resonance 



to higher frequencies, and to increase the number of bub- 
bles per unit volume to magnify their contribution. This 
explains the difference between the fitted and the actual 
size distribution in figure [5] inset. An important practi- 
cal application is that if one uses Foldy's model for fitting 
the experimental data, an incorrect bubble size distribu- 
tion will be obtained. Note that the same phenomenon 
was observed for bubbles in bread dough, with the same 
trend of Foldy's model predicting bubbles smaller than 
those determined from x-ray measurements [23j ] . 

Another possible source of error was the shear modu- 
lus, whose measurement was not very precise. For each 
frequency, we determined the complex shear modulus 
that gave the best fit to the experimental data. As shown 
by the triangles in Fig. [SJ a good fit was possible for the 
attenuation. For the velocity, the agreement was good for 
lower frequencies, but poor at higher frequencies. More- 
over, as shown in figure |51 the fitting leads to unrealistic 
values for \i' and n" . Again, one can qualitatively un- 
derstand these values: the real part of the fitted shear 
modulus is high in order to shift the peaks to higher fre- 
quencies, and the imaginary part is low (even zero) so 
that the oscillations of the bubbles are not damped by 
viscosity, and the scattering function / s is large. 

These results suggest that Foldy's model is inadequate 
for describing wave propagation in the PDMS sample. In 
the following we examine two possible explanations, and 
propose a model to account for these discrepancies. 
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FIG. 6: ( Color online) Shear moduli as functions of frequency. 
The lines correspond to the results for experimental [i! (solid 
line) and p!' (dashed line), as given in equations f!3\) . Open 
triangles show the fitted values for fj,' , solid ones for fi" . 



B. Homogeneity of the sample 

An important hypothesis made in all the effective 
medium models consists in assuming that the sample is 
homogeneous, i.e., any part of the sample can be con- 
sidered as equivalent to the others. The x-ray inspection 
of the sample reveals a homogeneous structure in the 
transverse directions, but a non-homogeneous one along 
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the thickness. Probably because of the upward buoy- 
ancy force exerted on the bubbles during the filling of the 
cell, the concentration of bubbles is varying, as shown in 
Fig. [7] the bottom part of the sample is free of bubbles, 
whereas the top part has a higher concentration. The 
bubble size distribution is also slightly different from one 
layer to the other. 

We investigated the influence of this heterogeneity of 
the sample on the acoustic propagation. Each layer i 
(1 < i < 50) of thickness d = 52 /xm was considered as an 
effective medium following Foldy's model with a concen- 
tration n'/d and a bubble size distribution as measured 
by the x-ray tomography. Following Brekhovskikh f2~ij ). 
we computed the transmission through this multi-layer 
system, and compared it to the experimental transmis- 
sion. Figure [5] shows the result: accounting for the inho- 
mogeneity of the bubble concentration with the 50 layer 
model does not change the prediction much, implying 
that the heterogeneity of the bubble concentration plays 
a negligible role in acoustic wave propagation through 
the medium. 




0.1 0.2 0.1 0.2 0^3 0.4 
ri (1/mm 2 ) R (mm) 

FIG. 7: (Color online) The layered structure of the sample is 
revealed by the x-ray tomography. Along the thickness h of the 
sample, we report the number of bubbles per unit area n' (left) 
and the radii distribution histogram (right) for each layer. 



C. Positional correlations 

Another approximation in Foldy's model is the as- 
sumption that the positions of the scatterers are uncorre- 
cted. Recent experiments in a 2D system of steel rods in 
water have shown that Keller's approach results in a very 
good correction to the model when correlations needed 
to be included [25|. In Keller's approach, the effective 
wave vector is related to the correlation function by: 

k 2 = k 2 + Aimfs 

- (4irnf s ) 2 J(l - g( r ))^^e ik ° r dr (15) 

where the correlations are taken into account by a radial 
function g(r) such that the local concentration of scat- 
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FIG. 8: (Color online) Amplitude (top) and phase (bottom) 
of the acoustic transmission through the sample. Because the 
data are the same as in figure^ the same discrepancy between 
Foldy's prediction for a homogenous sample (solid lines) and 
the experimental data (circles) is obtained. The prediction for 
a layered medium (dashed lines) does not bring better agree- 
ment. 



terers at a distance r from a scatterer is n(r) = nog(r) 7 
where no is the average concentration. For point-like 
scatterers, g = 1 and Foldy's approximation is recovered. 
For hard spheres, g = for r < 2R and 1 at larger dis- 
tances. In general, g can take more complicated forms, 
depending on the correlation mechanisms involved in the 
medium. 

From the x-ray tomography data, we were able to es- 
timate the function g for our sample. The general proce- 
dure for measuring such a function involves the following 
steps. First, one considers a sphere of radius rjy around 
one bubble and counts the number of bubbles in each 
spherical layer between r, and r,+i. Then, averaging 
over different central bubbles gives a statistical estimate 
of g(r). However, a limit of this technique is that the 
bubbles chosen as centres of the spheres cannot be at less 
than r^v from an edge of the sample. Thus, if one wants 
to investigate long-range correlations by taking high val- 
ues of r/v, the statistics are poor because the number of 
bubbles that are far enough from the borders is small. 
For our thin sample, this limitation was severe. To cir- 
cumvent this problem, we used 1/8 of each sphere, with 
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FIG. 9: (Color online) Function g as measured from the x- 
ray tomography data (circles). Error-bars correspond to one 
standard deviation. The solid line corresponds to the hard- 
sphere approximation for the bubble size distribution measured 
in the sample. 



an orientation chosen so that its volume did not cross 
any boundaries of the sample. Figure [9] shows the re- 
sult of these measurements up to a maximum distance of 
0.9 mm. It appears that the g function is very similar to 
the predictions of the hard-sphere approximation (solid 
line). Note that the oscillations in the g function that 
are predicted by the frequently used Perkus-Ycvick ap- 
proximation |25| are not expected to be significant here 
because the concentration of scatterers is low. Numer- 
ical simulations with random positions of scatterers in 
a box confirm that a simple hard-spheres law is a good 
approximation for a 2 % volume fraction. 

The predictions of Keller's model for the attenuation 
and velocity in the sample are shown in Fig. [TU]for three 
different g functions. When the hard sphere approxi- 
mation is used (i.e., the solid curve of Fig. [5] is taken), 
one obtains the dashed lines, which only give a slightly 
better agreement than the ISA. From the measured cor- 
relations, we also consider the high limit (i.e., the exper- 
imental points plus one standard deviation) and the low 
limit (i.e., the experimental points minus one standard 
deviation) for the g function. They give different cor- 
rections: while the former gives a worse agreement, the 
latter improves it. One can understand this difference 
qualitatively by re-writing equation (|15[) as 



f 



47m/ s 

k 2 
o 



1 — 47m / s J (1 — g)rdr 



(16) 



where the reasonable approximations kr <C 1 and k^r <C 
1 have been made. It is important to note that, in the 
frequency range over which the discrepancy is high, the 
bubbles are responding in phase opposition to the pres- 
sure field, i.e. Re(/ S ) < 0. So, according to equation (fT6l) . 
the correlations increase the effect of the bubbles on k if 
g < 1, whereas they decrease it if g > 1 . Note that a sim- 
ple criterion can be deduced from Eq. (|16|) to estimate 
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FIG. 10: (Color online) Investigation of the effect of corre- 
lations. The experimental attenuation and velocity (circles) 
and the ISA prediction (solid thin lines) are shown again, as 
in figure^ Keller's approach (Eq. is plotted for the 

hard-sphere correlation (dashed lines), for the high estima- 
tion of g measured from the x-ray data (dotted lines), and for 
the low estimation of g (thick solid lines). 



the importance of the corrections due to correlation. If 
g is roughly described as a step function going from to 
1 at distance r c , the term with the integral in Eq. (TH 
which we will denote as /3, reduces to 



2 R \RJ 



(17) 



Thus, the magnitude of the correction depends on the 
ratio of r c to the radius of the bubbles R. Long-range 
correlations are thus expected to have a stronger effect. 

It appears that Keller's approach, with a reasonable 
estimate of the g function, introduces a correction that 
goes in the right direction for better agreement. But 
the agreement is still unsatisfactory. Moreover, Keller's 
expression applies to a monodispcrsc medium, whereas 
our bubbly PDMS is a polydisperse sample. So far, 
we have assumed that the polydispersity could be taken 
into account by changing nf s into J n(R)dRf s (R) in the 
equations. However, as we show in appendix a self- 
consistent approach provides a more general expression 
(see equation (|A7[) ). in which the correlations depend on 
the size of the scatterer: g(r) becomes g(r, R). To exam- 
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FIG. 11: (Color online) Function g(r,R) as measured from 
the x-ray tomography data ( circles ) on three populations of 
bubble radii. Error-bars correspond to one standard deviation. 
The solid line corresponds to the hard-sphere approximation 
for the bubble size distribution measured in the sample. 
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ine this polydispersity effect, we considered 3 different 
populations of "central" bubbles: 0-140 /zm, 140-180 /zm, 
and 180-600 fim. For each population, we were able to 
average over approximately 1000 bubbles, ensuring good 
statistics. Figure QT] shows the g functions obtained for 
these 3 populations. For the small bubbles, there is a 
clear maximum of the g function. This is an indication 
of clusters. A close inspection of the picture in figure [H 
indeed gives the impression that the small bubbles are 
clustered. For the average size bubbles, the g function 
is similar to the hard-sphere approximation. For the big 
bubbles, we find evidence of depletion: there is a lower 
probability of finding a bubble close to a big one. Again, 
note that the photo in Fig. [5] shows isolated big bubbles. 

Figure [12] shows the predictions of our self-consistent 
approach (SCA) when 3 different possibilities are consid- 
ered for g(r,R): hard sphere approximation (solid lines 
in Fig. QT]), an d the measured values of g including the 
low and the high limits (bottom and top of the error 
bars in Fig. [TT]) . It appears that the effects observed in 
Fig. [TU] arc amplified when the polydispcrse correlation 
is taken into account. Interestingly, even for the hard 
sphere approximation, slightly better agreement is found 
(see table ITO] ). This can be explained by the following 
argument: large bubbles scatter more than small bubbles 
(see equation ((3])), and their average correlation distance 
r c is also larger. 

When the range of measured values of g is considered, 
the agreement is better or worse, depending on whether 
the high or the low limit of the possible g function is 
considered. It follows that, given the precision of our 
correlation measurements, we cannot conclude whether 



FIG. 12: (Color online) The effect of polydisperse corre- 
lations as predicted by the self-consistent approach (SCA). 
Three g(R,r) functions are considered (see figure hard- 
sphere correlation (dashed lines), high and low estimates of g 
measured from the x-ray data ( dotted and thick lines, respec- 
tively). 



TABLE III: A simple criterion for estimating how well the 
predictions of the different models agree with the experimen- 
tal data: the relative mean squared differences between the- 
oretical and experimental effective wavenumbers (average of 
the imaginary and real parts) over the 30-300 kHz frequency 
range. Different g functions are considered: from the hard- 
sphere approximation (hs), and from the measured correlation 
(see Fig. Hit with the lower (xp-) middle (xp) and upper (xp+) 
estimates. 



Model 


ISA 


Keller 


SCA 


9 


/ 


hs xp - xp xp + 


hs xp - xp xp + 


(%) 


14.8 


11.3 7.3 12.4 17.7 


10.1 4.3 10.8 19.6 



or not the polydisperse correlations are sufficient to fully 
explain the discrepancy. Nonetheless, we find that plau- 
sible values for g are able to give good agreement. A 
more precise determination of the correlations would be 
necessary for a definitive conclusion to be drawn. 
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V. CONCLUSIONS 

Wc have investigated wave propagation through a com- 
plex medium in which low-frequency scattering reso- 
nances lead to strong dispersion, with large peaks in the 
velocity and attenuation. By choosing to work with a 
relatively simple acoustic system - bubbles in the elastic 
medium PDMS - we have shown how these signatures of 
strong resonant scattering are influenced by the coupling 
between the scattcrers. This leads to the breakdown of 
the independent scatterer approximation (ISA), which is 
commonly used to interpret experimental data. A cru- 
cial step in investigating the failure of the ISA in this 
system has been our ability to characterize the physical 
properties of bubbly PDMS very carefully, including in- 
dependent measurements of the elastic properties of the 
matrix as well as the concentration and size distribution 
of the scatterers. Remarkably, we find that this failure 
of the ISA occurs even though the concentration of scat- 
terers is low, only 2% in our case. It is worth noting 
one important consequence: if the ISA is used to esti- 
mate the size distribution of the scatterers, significantly 
smaller sizes are found than by direct imaging methods. 

In this paper, we have proposed that this discrepancy 
can be explained by accounting for the role of correlations 
in the scatterers' positions, which we were able to probe 
directly using x-ray tomography. Following Keller's ap- 
proach, we have shown that these correlations brought a 
non-negligible correction to the ISA, and that the agree- 
ment with the experimental data was improved, although 
not perfect. We have proposed an alternative approach, 
based on a self-consistent argument, with which we were 
able to take into account the effect of polydispersity on 
the correlations: bubbles of different radii are not cor- 
related with their neighbours in the same way. Again, 
x-ray tomography allowed us to investigate this polydis- 
perse correlation, which was found to be different from 
the monodispersc estimate: small bubbles are more likely 
to be close to other bubbles, whereas big bubbles are of- 
ten isolated. We have shown that, within the possible 
polydisperse correlations that are consistent with our x- 
ray measurements, the SCA was able to provide a satis- 
factory explanation for the observed shift to higher fre- 
quencies in the acoustic resonances in the sample. It 
should also be noted that our model could be refined to 
account for more complex situations. In particular, when 
the correlations are long-range in a thin sample, one can 
expect the boundaries to play a role in the correlation, 
an effect we did not incorporate. Also, we looked at cor- 
relation effects or inhomogencity effects, but we have not 
considered the case in which both effects are coupled. 

An interesting question is the following: why are cor- 
relations important in the bubbly medium studied here, 
whereas they have not been detected before for bubbles in 
water [2|| or bubbles in a yield stress fluid Q? According 
to equation (jTTJ) , the magnitude of the correction due to 
the positional correlations is proportional to the volume 
fraction of scatterers 4>, which was generally lower in pre- 



vious studies (5 x 1CP 4 to 10~ 2 at most), so that j3 would 
be expected to be low. However, it is interesting to note 
that some previous studies did report deviations between 
Foldy's prediction and the experimental data that might 
be interpreted as a manifestation of correlation effects. 
In the historical data by Silberman [27|, for instance, the 
measured attenuation was much larger than Foldy's es- 
timate. Feuillade proposed that the discrepancy could 
be resolved by introducing an arbitrary fitting param- 
eter [28| that one can interpret as a correlation length 
(see the appendix). Wilson et al. also observed devi- 
ations from Foldy's estimate, as their attenuation peak 
was at a lower frequency than what was expected from 
the size analysis (see Fig. 10 in !2G|). They attributed 
this discrepancy to the uncertainty in their measurement 
of the bubble size distribution, but one can imagine that 
correlation effects were responsible for this shifting of the 
attenuation peak. Indeed, large depletion zones around 
bubbles have been observed in bubbly liquids [29[ , which 
would lead to positional correlations. In our sample, the 
exact mechanism leading to the positional correlations 
has not been identified unambiguously. Ostwald ripen- 
ing might induce the disappearance of small bubbles close 
to bigger bubbles, hence leading to the isolated big bub- 
bles we observed. Another possible scenario could involve 
the rotation of the cell when the sample was prepared. 
If one considers that only the buoyancy and the Stokes 
drag force are relevant (inertial forces were negligible in 
our set-up), a bubble is expected to move along a cir- 
cle whose radius is inversely proportional to the rotation 
frequency, and proportional to the square of the bub- 
ble radius. The movement of the bubbles may induce 
hydrodynamic interactions leading to positional correla- 
tions. This hypothesis is supported by the observation 
of strong clustering when the cell was rotated at a lower 
speed (0.7 rpm), i.e., when the radii of the circles followed 
by the bubbles were larger. 
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Appendix A: Self-consistent approach for multiple 
scattering of waves 

We propose a self-consistent approach for calculating 
the dispersion relation of a random multiply scattering 
medium, taking into account the positional correlations 
of the scatterers. Let us consider an infinite medium with 
n inclusions per unit volume. We limit ourselves to low 
frequencies such that the wavelength is large compared 
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to the size of the scattering inclusions, which means we 
can consider that the waves scattered by the inclusions 
are spherical: when scatterer i is excited by a wave with 
amplitude pi, it generates a wave Pifl exp(ifco?")/r a t dis- 
tance r, where /* is its scattering function and ko the 
wave number in the pure medium. The total field ex- 
perienced by inclusion i, due to the other scatterers, is 
given by the self-consistent relation 



» • 




X 



(Al) 



where is the distance between scatterers i and j. In 
the case of a finite number of scatterers, Eq. (|Alj) can be 
solved and all the pi exactly computed. For an infinite 
number of inclusions, however, taking into account the 
full multiple scattering is a difficult task and the pi are 
not easily accessible. In particular, the 1/r range of the 
interaction means that no pair can be neglected a 
priori in Eq. (|Aip . 

Let us further simplify the problem by assuming that 
all the scatterers arc identical (/* = / s ) and look for 
modes that take the form of plane waves propagating in 
the x direction: pi = Pexp(ifcxi), where k is the effective 
wave number we want to determine, and P the amplitude 
of the mode. Then, Eq. (|A1[) reduces, for scatterer i = 
arbitrarily chosen as the central one located at r = 0, to 



1 



#0 



(A2) 



where rj is the distance to scatterer j, and 9- } the angle 
between the direction of propagation and the position of 
scatterer j (see Fig. fT5|) . 

One can approximate the discrete equation by a con- 
tinuous equation, taking into account the positional cor- 
relations with the pair correlation function g(r): 



f s J d 3 r x ng(r) 



^(krcosG+kor) 



ng(r) sin(/cr)e ° r dr. 



(A3a) 
(A3b) 



Invoking a small imaginary part of k to preserve the 
convergence of the integral, one can then calculate the 
effective wave number: 



k 2 



» n 



Aimf s 



1 + Annf s J(l- g(r))^^c ik ° r dr ' 



(A4) 



Note that equation (|A4[) reduces to Foldy's equation if 
the denominator is approximated by 1, and to Keller's 
equation if it is expanded to second order in nf s . Interest- 
ingly, it has some similarities with the equation proposed 
by Feuillade, which was found to give excellent agreement 



FIG. 13: (Color online) An infinite 3D medium with spherical 
scatterers is considered. The position of scatterer is arbi- 
trarily chosen as the centre of the coordinates and we look for 
modes of the system consisting of plane waves propagating in 
the x direction. The dispersion relation is obtained by calcu- 
lating the total field experienced by scatterer due to the other 
scatterers. 



with experimental data for bubbly liquids (see equation 
(30) i n I3CI1. Feuillade's model has been subject to criti- 
cisms [13ll3l|. in particular because, for unclear reasons, 
it assumes that the interactions between the scatterers 
are confined within a finite range. It might be that this 
finite range is related to correlations. Note that the ap- 
proach we develop here is similar to Feuillade's. The 
differences are that (1) we abandon the assumption that 
all the scatterers are in phase, which is true only locally, 
and (2) we include the positional correlations. 

For a polydisperse assembly of scatterers, the same 
self-consistent scheme can be used. In this case, the total 
field pi depends a priori on the radius of the bubble, 
Pi = P(Ri)exp(ikxi), meaning that equation (|A3a[) must 
be generalized into 



P(R) = f s (R) J n{R')AR'P{R') 



g(r,R,R') 



^ i ( k r c o s + k o r ) 



d 3 r, 



(A5) 



where g(r 7 R, R') gives the probability of finding a bubble 
with radius R 1 at a distance r from a central bubble with 
radius R. In the general case, equation (|A5|) cannot be 
directly simplified into a dispersion relation. However, 
if we assume that the correlation does not depend on R' 
(g(r, R 7 R') = g(r, R)), a further integration over R yields 

n(R)f s (R)AR J g(r,R) d 3 r(A6) 

from which we obtain the dispersion relation 



2 2 jAnn(R)f a (R)dR 



(A7) 



with 



A = / 4nn(R)f s (R)dR / (1 - g(r,R)) 



sin kr 



ik r 



k 



dr. 
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Note that, as in the ISA, polydispersity modifies the 
equation by changing nf s into J n(R)dRf s (R). But poly- 
dispersity also brings another non-trivial modification: 
the correlation can depend on the radius of the cen- 



tral bubble considered, which makes different correlations 
possible for different bubble sizes, as in our bubbly sam- 
ple. 
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